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One of the famous results of network science states that networks with heterogeneous connectivity are more 
susceptible to epidemic spreading than their more homogeneous counterparts. In particular, in networks of iden¬ 
tical nodes it has been shown that network heterogeneity, i.e. a broad degree distribution, can lower the epidemic 
threshold at which epidemics can invade the system. Network heterogeneity can thus allow diseases with lower 
transmission probabilities to persist and spread. However, it has been pointed out that networks in which the 
properties of nodes are intrinsically heterogeneous can be very resilient to disease spreading. Heterogeneity in 
structure can enhance or diminish the resilience of networks with heterogeneous nodes, depending on the cor¬ 
relations between the topological and intrinsic properties. Here, we consider a plausible scenario where people 
have intrinsic differences in susceptibility and adapt their social network structure to the presence of the disease. 

We show that the resilience of networks with heterogeneous connectivity can surpass those of networks with 
homogeneous connectivity. For epidemiology, this implies that network heterogeneity should not be studied in 
isolation, it is instead the heterogeneity of infection risk that determines the likelihood of outbreaks. 


In the exploration of complex systems, epidemiology plays an important role as a source for toy models and case studies, but 
also an area where a real world impact can be made 11 - 6]. It has been pointed out that new diseases have emerged whenever 
environmental change brought humans in contact with new pathogen or disease vectors, i.e. animal hosts of a given disease 0]. 
The past decades have brought rapid environmental change, a growing world population, and increasing long-range connectivity 
in relevant networks, due to human travel and livestock transports 0. Together with decreasing vaccination levels and misuse of 
antibiotics, this has led to both emergence of new diseases and the return of old ones, sometimes in the form of highly resistant 
strains. In the future antibiotics, vaccinations, and quarantine are bound to remain first line of defence against these threats. 
However, insights from physics that can improve the efficiency of these measures, even if only by a small measure, have the 
potential to save many lives and relieve the economic burden created by the disease. 

Many current studies seek to determine the so-called epidemic threshold, the critical level of infectivity that a pathogen needs 
to surpass to spread and cause large outbreaks This threshold depends on many factors including the structure of the 

underlying network of contacts, the heterogeneity in the hos t po pulation, and behavioral responses to the disease. Among these, 
the effect of network structure is perhaps best understood 111 ll4l4fl . It can be shown that the ability of a disease to spread is 
generally related to the leading eigenvalue of the networks adjacency or non-backtracking matrices js| 13-18]. Factors that 
increase this eigenvalue, lead to lower epidemic thresholds. Two well-known factors that facilitate the spreading of diseases are 
high network connectivity and network heterogeneity. Here, network heterogenity can be defined as the second moment of the 
degree distribution, the probability distribution of the number of links on a randomly chosen node. In extreme cases, for example 
in scale-free networks, the epidemic threshold may vanish in the thermodynamic limit such that diseases with arbitrarily low 
infectivity can still spread lill[ |19fl . 

With respect to individual heterogeneity in the host population, some questions remain open. Generally intra-individual and 
link heterogeneity 0, [2lll , such as different levels of resistance, times of contact, differences in infectitivity or recovery rates, 
reduce the size and risk of epidemics 12214271. For instance, it was shown analytically and numerically that epidemics are most 
likely if infectivity is homogeneous and least likely if the variance of infectivity is maximized 1221 12311 . Comparable results 
were found in lattice models and in a biological experiment |24l| . However, heterogeneous susceptibility can make networks 
more vulnerable to the spread of diseases if the correlation between a node’s degree and susceptibility are positive 12711 . This is 
intuitive as the intra-individual heterogeniety, in this case, amplifies the negative effect of network heterogeniety. Finally, it was 
reported that heterogeneous susceptibility can in some scenarios cause a secondary epidemic after a primary outbreak 1 2cl 281. 

The studies referenced above focussed on epidemics on static networks. Another active line or research concerns the effect of 
dynamic, adaptive 12^- 32 ] or temporal networks 133, 34]. While these types of networks are closely related, dynamical networks 
emphasize the overall statistical effect of changing connectivity 12911301, adaptive networks emphasize the dynamical response 
of network structure to the disease state[351136], and temporal networks focus on the impact of the specific timings of events 
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An area that is so far poorly understood is how the behavioral response of individuals to epidemics reshapes the network. 
Beside institutionalised responses, such as mandatory vaccinations and quarantine, humans react to the outbreak of a major 
epidemic in a variety of ways, for instance by increasing hygiene, using protective measures such as face masks, reducing 
contact with infected individuals l39ll - or avoiding contacts with other humans when infected. The common element in these 
responses is a reduction in the frequency of contacts between infectious and susceptible individuals. This limits the disease 
propagation by decreasing the effective network connectivity. However, the ultimate effect of behavioral responses cannot be 
understood as a static reduction of connectivity. If individuals respond to the epidemic state of other individuals by altering their 
interactions then a complex dynamical feedback loop is formed in which the state of nodes affects the evolution of the topology, 
while the topology governs the dynamics of the nodes’ state. Thus an adaptive network is created. 

In the context of epidemics it has been shown in a simple model [40Q that trying to avoid contact with infected is highly 
effective against small outbreaks but less effective against an established epidemic. Instead of a single epidemic threshold, such 
models possess an invasion threshold for new diseases and a (lower) persistence threshold for established epidemics. Thus, 
a parameter region is formed where a disease cannot spread when it is newly introduced, but can persist when it is already 
established. Physically speaking, a hysteresis loop is formed and the percolation transition at the onset of the epidemic becomes 
discontinuous. Under certain conditions the behavioural feedback loop can also lead to the emergence of epidemic cycles 


[35 


The analysis of wide variety of related models showed that these observations are robust over a wide class of models £10- 
451. Furthermore, studies showed that behavioral response can increase the impact of targeted vaccinations |46||. and that the 


timing of interventions can be more important than in static networks [47]. 

In this paper, we study the combined effect of heterogenity in intra-individual parameters and the behavioral response. We 
show that, starting from a well-mixed network, a heterogeneous connectivity is formed. It is known that heterogeneous networks 
of heterogeneous nodes can be very resistant to disease outbreaks, if certain correlations are present|27]. Here we show that these 
correlations naturally arise in the adaptive network, and that the resulting network configuration is generally significantly more 
resistant to outbreaks than a network with homogeneous topology. Our analysis suggests that the decisive property governing 
disease invasion is not network heterogeneity but the heterogeneity of the effective disease risk of agents. 


Results 


Heterogeneous adaptive SIS model. We consider a network of N agents connected by K bidirectional links. We distinguish 
two types of agents, which we denote as type A and type B, which differ by their resistance to the disease. The type is an internal 
property of the agent that does not change. Hence, the proportion of agents of type A, p a , and the proportion of agents of type 
B, pb = 1 — Pa, are constant in time. For conciseness we denote the total number of agents in a given state i G {a, b} by 
Ni = Npi. Moreover, nodes carry an additional internal variable that indicates their epidemic state. We consider a variant of the 
susceptible-infected-susceptible (SIS) model of epidemic diseases (48], such that a given node is either susceptible to the disease, 
state S, or infected (and infectious), state I. The proportion of nodes in the S and I state is denoted as [5] and [7] = 1 — [5], 
respectively. 

The network is initialized as an Erdos-Renyl random graph, G(N, K). These networks have a narrow, Poissonian degree 
distribution, such that the network connectivity is homogeneous in the initial state. Each node is randomly assigned a type and 
epidemic state such that desired values of p, and the desired initial values of [,Sj and [7] are realized. Time evolution of the 
network is then driven by three processes, namely the a) recovery of infected nodes, b) contact avoidance, and c) contagion. 
These are implemented as follows: a) Infected nodes recover at rate p, returning to the susceptible state, b) A given link, 
connecting a susceptible agent to an infected agent is rewired at rate w. In an rewiring event the original link is cut and a new 
link between the susceptible node and a randomly chosen other susceptible node is created, c) For every link connecting a 
susceptible to an infected node, the disease is transmitted along the link at a rate /3ipi that is dependent on the type i of the 
susceptible node. 

In the following we assume that ip a > ipb such that nodes of type A are more susceptible to the disease than nodes of type B. 
Our mathematical results hold for parameters in arbitrary units, but the rates can be thought of having the dimension of nodes- 
per-time (recovery) and links-per-time (contagion, rewiring) respectively. Throughout the paper we balance the parameters ip a 
and ip b such thatpa^a + Pbtpb = (tp)- The variance of susceptibility is crj = ((-0) — ipb)(tpa — (0)). 

Most of our results below are found by analytical calculation or continuation of solution branches and are thus non-simulative 
in nature. However, we compare these results to large agent-based simulations. In these simulations we use an event-driven 
(Gillespie) algorithm to simulate the stochastic process described above. In simulation we use an Erdos-Renyi network with 
N = 10 5 nodes and K = 10 6 links with recovery rate p = 0.002, rewiring rate lo = 0.2 and average susceptibility (ip) = 0.5, 
unless noted otherwise. We vary ip a as a proxy for heterogeneity. For every choice of ip a and iph, the parameter p a is set such 
that (ip) = 0.5. All parameters used in simulation runs are stated in the caption of the respective figure. 
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FIG. 1: (Color online) Bifurcation diagram of the adaptive heterogeneous SIS model. Show is the stationary level of disease prevalence 
I* as a function of infectivity /3. When the infectivitiy is decreased the endemic state vanishes in a saddle node bifurcation (‘Persistence 
threshold’). The disease free state can be invaded by the epidemic, if the infectivity surpasses a point where a transcritical bifurcation occurs 
(‘Invasion threshold’). Agent-based simulation (circles) and equation-based continuation (lines) provide consistent results on the persistence 
threshold, but predict different invasion thresholds. This discrepancy appears due to a projection effect, because the state where the disease 
is extinct is not uniquely defined (see text). In addition to stable solution branches (solid) the continuation also reveals an unstable solution 
branch (dotted). Parameters: i /) a = 0.65, t/)b = 0.05, p a = 0.75, io = 0.2, /x = 0.002, N = 10 5 , K = 10 6 . 


Hysteresis in the heterogeneous model. To gain some basic intuition, let us first investigate the system by explicit agent- 
based simulation of the network model |49j|. For this purpose we evolve the system in time, according to the stochastic rules, 
until a stationary level of epidemic prevalence is approached. Repeating this procedure for different values of infectivity 3 
reveals the diagram shown in Fig. Q] The Figure is qualitatively similar to results from the homogeneous adaptive SIS model: 
Epidemics starting from a small proportion of initially infected agents go extinct deterministically if the infectivity is below a 
certain threshold, /3i nv , which we identify as the invasion threshold. By contrast, established epidemics, i.e. simulation runs 
starting with a higher initial proportion of infected, can persist if the infectivity surpasses a different persistence threshold /3 per - 
The persistence threshold is lower than the invasion threshold, such that a bistable region is created. In this region an epidemic 
that enters the system at low density goes extinct, but an established large epidemic, that perhaps entered the system earlier when 
parameters were different, can persist. The bistable region constitutes a hysteresis loop: If we slowly increase the infectivity 
the extinct state is stable until where a jump to the endemic state occurs. If we lower /3 again, the system persists in the 
endemic state up to /3 per , where it collapses back down to the extinct state. 

In the following we explore the effect of heterogeneity on the thresholds /3; nv and /? per . To gain analytical insights, we consider 
a moment expansion of the system |50]. We use symbols of the form [X. u \ and [X U Y V ] with A', Y G {I, S} and u, v G {a, b} to 
respectively denote the proportion of agents and per capita density of links between agents of a given type. For instance [/ a ] is 
the proportion of agents that are infected and of type A, and [5 a /b] is the per capita density of links between susceptible agents 
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of type A and infected agents of type B. All of these variables are normalized with respect to the total number of nodes N. Given 
the number of infected nodes of a given type we can thus find the number of susceptible nodes by using the conservation law 

[I u \ + [S u \=Pu. 

The time evolution of the proportion of nodes that are infected and of type A and B can be respectively written as 

^-[4] = -m[4] +0ip a ^[5 a J„], (1) 


±[i h ] = - fl [i b ]+^ h j2[s h i v }. 

V 

For the link densities, using a pair-approximation leads to equations of the form 


d[S a S a } 

ctt 


=ll[S a I a \ - 2^a( 


[5 a 5 a ][S a / a ] 

[5a] 


[5 a 5 a ] [5 a / b ] 
[5a] 


) + 


w[5 a 


[5a] + [5b] 


([5 a /a] + [5 a / b ]), 


( 2 ) 


(3) 


where the terms on the right hand side describe the impact of the different processes on the motif considered, [jS a jS a ] in this 
example. For instance the first term corresponds to the creation of S' a -S' a -links due to recovery of the infected node in S a -I a - 
links. In total the / a -nodes recover at the rate /i[/ a ]. Every such recovery event creates an expected number of S'a-S'a-links that 
is identical to the average number of / a -S a -links anchored on an / a -node, which is [/ a S a ]/[/ a ]. In summary, the change in the 
density of S a -S a -links due to recovery of / a -nodes is /i[/ a ][4S a ]/[J a ] = /i[/ a S a ], which explains the first term in Eq. (0). 

In total the moment expansion yields a system of 11 ordinary differential equations. For conciseness we show the remaining 
equations in the Methods. 

We solve the moment equations by numerical continuation of solution branches using AUTO [|5ji]. This reveals branches of 
stable and unstable steady states (Fig. 1). As in homogeneous adaptive SIS model the limits of the hysteresis loop are marked by 
a fold bifurcation and a transcritical bifurcation point. In the fold bifurcation the endemic steady state collides with an unstable 
saddle and the two states annihilate. In the transcritical bifurcation the saddle state intersects the healthy steady state, which 
causes the healthy state to become unstable. The value of [3 in the fold bifurcation point thus marks the persistence threshold 
/3 per and the critical value in the transcritical bifurcation point marks the invasion threshold /3; nv . 

The comparison between the continuation results and agent-based simulation, in Fig.0 shows that both methods are in good 
agreement regarding the location of the solution branches. Also the values for the persistence threshold agree. However, the 
continuation predicts a much higher value of the invasion threshold than the agent-based simulation. 

To understand the discrepancy in the thresholds let us again consider the plot in more detail. The diagram shows a projection 
of the full 11-dimensional space spanned by the moment equations. In general bifurcation diagrams of this type identify the 
steady states uniquely. However, this is not the case when the epidemic is extinct. If there are no infected left, then the dynamics 
freezes independently of the connectivity of the nodes. Thus, as long as [J a ] = [J b ] =0 the state under consideration is stationary 
regardless of the values of [<S' a S' a ] and [St, St,]. Hence the zero line in the bifurcation diagram is really a 2-dimensional plane of 
absorbing steady states. 

We can now resolve the discrepancy between the continuation and the simulation results. Continuation of the unstable solution 
branch leads to a specific point on the plane of extinct states. This landing point is uniquely defined and marks the invasion 
threshold of the network with the corresponding values of [S a S a ] and [S b S b ]. However, these values are not identical to those 
considered in the numerical simulation. 

We argue that both of the invasion thresholds have significance. The simulation is valid for the well-mixed initial system, 
where the network structure has not yet adapted to the presence of the disease. Thus this threshold is relevant in case of the 
arrival of a new disease. By contrast the threshold found by continuation corresponds to a case where the network structure has 
adapted to the disease, for instance due to repeated exposure to same or similar pathogens. In the following, we therefore refer 
to the two thresholds as the initial and adapted invasion threshold, respectively. 

Invasion thresholds and heterogeneity. We emphasize that the observed discrepeancy between the initial and the adapted 
invasion threshold could not appear in networks of identical nodes. For identical nodes the extinct state is unique on the level 
of the pair approximation ([/] = [II] = [SI] = 0), and thus both thresholds must coincide. The results in Fig. [2] show that is 
indeed the case, while different thresholds are observed in all networks with heterogeneous nodes. 

We note that the adapted invasion threshold is always higher than the initial threshold. Thus adaptation increases the robustness 
of the system to disease invasion. Let us therefore explore the adaptation in more detail. The adaptation is driven by the rewiring 
which means that nodes that are frequently infected in average lose links, while nodes that are rarely infected gain links. On the 
population level this means that the average degree of nodes of type A, fc a , and the average degree of nodes of type B, fc b change 
dynamically in response to the average proportion of nodes of each type that are infected. Before trying to compute the ratio 
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FIG. 2: (Color online) Comparison of thresholds. The plot shows a very good agreement agreement between equation-based continuation 
(lines) and agent-based simulations (symbols) for the persistence thresholds /3 per (box, dashed). However, a notable difference exists for the 
invasion thresholds /3i nv (circle, dotted). Parameters: ipb = 0.05, ui = 0.2, fi = 0.002, N = 10 5 , K = 10 6 . 


fcb/fca let us first point out that a lower bound is kh/k :l = 1, in this case the degrees are equal, and thus nodes of type A would 
be infected more frequently, due to their higher susceptibility. Hence k H would decrease and the ratio k\,/k H would increase. An 
upper bound is provided by fcb/£: a = V’a/'^b- In this case, i/Ja.k a = ipbkb implies the that both types of nodes get infected at 
equal rate. Because the degree of nodes of type B is higher than the degree of nodes of type A, an infected node of type B will 
lose links more rapidly and hence the ratio k\ t /k- rl will decrease. 

The numerical value of the degree ration k\Jk rL is shown in Fig. [3] To gain also an analytical understanding we resort to a 
description of the system that is coarser-grained than the full moment expansion. First, note that, on a population level fc, ; , the 
mean degree of nodes of type i £ {a, b}, obeys a differential equation of the form 



-ki{Ii]u + [Sjju, 


(4) 


where the two terms denote rewiring losses and gains, and auxialliary variables u and v have been defined to contain all other 
factors which do not depend on the index i. In the steady state we find 


Thus, when we compute the ratio 


ki[Ii] = [Si]-. 

u 

fc b [4] [5 b ] 
M4] [Sa] 


(5) 


( 6 ) 
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FIG. 3: (Color online) Network heterogeneity in the adapted state, indicated by the degree ratio . The dependence of the degree 
ratio in agent-based simulations (black circles) closely follows the prediction from integration of the ODE model (solid red line), a very good 
approximation is also provided by the relationship kh/k a = yVv/V’b (blue dashed line), whereas the naive expectation kb/k a = ipa./iph 
(pink dotted line) overestimates the network heterogeneity significantly. Contrary to expectations the networks following the naive solution 
(the most heterogeneous case), would be maximally stable against disease invasion. Parameters: tpb = 0.05, c o = 0.2, q, = 0.002, N = 10 5 , 
K = 10 6 . Inset: the magnification of the superimposed part of the main figure. 


the factors u and v vanish. Using a mean field approximation, the epidemic state variables follow equations of the form 

^ [Si\ = -qkiil>i[S 2 \ + r[Ii] (7) 

where the terms capture the effects of contagion and recovery, respectively, and again auxiliary variables q and r have been 
defined that contain all other factors that don’t explicitly depend on i. We use the same trick as before and consider the steady 
state, where 


qk^ilSi] = r[Ii] 

and hence 

kbtpbjSb] _ |7b] 
fcaV’al-S'a] [4] ’ 

Substituting this result into Eq. d6j we find 

fcb 2 ^b [5 b ] _ [5 b ] 
fcaVa [5a] [5 a ] 


( 8 ) 

(9) 


( 10 ) 
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FIG. 4: (Color online) Comparison of thresholds in self-organized networks. The plot shows a very good agreement between equation-based 
continuation (lines) and agent-based simulations started in an artificially created adapted state (symbols) for both the invasion thresholds /3i nv 
(circle, dotted) and the persistence thresholds /3 per (box, dashed). See Fig. 2 for comparison. Parameters: iph = 0.05, uj = 0.2, fi = 0.002, 
N = 10 s , K= 10 6 . 


and hence 


fcb 

K 



(ll) 


The result of this mean field argument is in good agreement with numerical results (Fig. 3). It implies that the rewiring mecha¬ 
nism considered drives the system to a state where the less susceptible nodes (type B) have a higher mean degree, which partly, 
but not fully, compensates for their lower susceptibility. Thus in the adapted network the nodes of type B get infected more often 
than they would in a network with homogeneous degree distribution, but still less often than the nodes of type A. 

To verify that the self-organization of the link distribution explains the observed discrepancy between the initial and the 
adapted invasion threshold, we turn to the agent-based simulation again. However, in this case we start the simulation from an 
artificially created adapted state. To initialize this state we simulate the system with the same set of initial parameters until the 
system reaches the stationary state. Then we retain the self-organized link pattern, but reassign all epidemic states, such that all 
agents are susceptible except for 20 initially infected. Then we simulate the system again until it either reaches the endemic state 
or the epidemic goes extinct. We locate the epidemic threshold by running a series of such simulations and find the point where 
the probability to reach the disease-free state becomes zero. The epidemic threshold that is thus found coincides with the result 
from continuation of the equation-based model (Fig. 0. 

Higher thresholds in heterogeneous networks. Let us emphasize that the adapted network has a more heterogeneous degree 
distribution (Fig. 0 but also a higher epidemic threshold. Thus comparing the adapted and the maximally random state, the 
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more heterogeneous network is actually more robust to the invasion of the pathogen. This appears counter-intuitive in the light 
of many recent work in network epidemiology, which showed that more heterogeneous contact networks are generally more 
susceptible to disease invasion. However, it was already shown in (27:1] that heterogeneous networks become less susceptible if 
there is a correlation such that the highly connected nodes have lower susceptibility to the disease. The present results show that 
a simple but plausible local adaptation rule can drive this process so far that the network with heterogeneous degree distribution 
is more robust against disease invasion than an Erdos-Renyi random graph of the same mean degree. 

Link heterogeneity is clearly a double-edged sword. On the one hand it is intuitive that some amount of heterogeneity in the 
degree is advantageous if it means that more links lead to nodes with low susceptibility. On the other hand, epidemic theory 
suggests thatTn the limit of very heterogeneous networks, the robustness of the network should decline because of the high excess 
degreeJHl 52|. This suggests that there should be an optimum level of heterogeneity, where the epidemic invasion threshold is 
highest. 

We can understand the interplay between the two effects by a link-centric percolation argument. We consider the limit of low 
disease prevalence and focus on the active links, i.e. links connecting an infected node to a susceptible node. Given a single 
focal active link we can estimate the expected number of secondary active that is created by transmission along the focal link by 


_ Pafc^a/3 Pb^bP 

Z/n — -=-1-=- 

pk pk 


( 12 ) 


where k = p a k a +pbkb is the constant mean degree of the system. Defining q = kb/k a and substituting k a = k / (p a + Pbq) and 
kb = k/{p a q~ l +Pb) we can express the link reproductive number Zq as a function of the degree ratio q. Although the resulting 
expression for Zq is relatively complex, we can find it’s minimum, i.e. the most robust point, by differentiating, which yields 


kb _ 0a 
k a t/>b 


(13) 


This coincides with the upper bound for the degree ratio, or, in other words, the point where the nodes of types A and B are 
infected at an identical rate. Therefore increasing the degree heterogeneity by increasing the degree ratio is advantageous to the 
point where the more resistant agents become infected more often than their less resistant counterparts. Thus it appears that the 
decisive characteristic that determines its robustness to disease invasion is not the heterogeneity of the the network structure, but 
the heterogeneity of the disease risk to which the agents are exposed. 

The independence of the result above from other parameters suggests that it is true in a wider class of systems, but this intuition 
will need to be validated in further investigations. For the adaptive system this result means that the network always operates in 
the regime where a higher degree ratio and therefore more heterogeneity has a stabilizing effect. We could in principle construct 
a system that self-organizes to the optimal point, by replacing the per-link rewiring rate by a rewiring rate per infected node. 
However this variant of the model is beyond the scope of the present paper. 


Discussion 

In this paper we studied an adaptive heterogeneous SIS model numerically and analytically. The analysis revealed that 
heterogeneity in the intrinsic parameters of the nodes induces heterogeneity in the connectivity: Over time more resistant nodes 
gain more links until a steady state is reached, in which nodes with higher resistance are still less likely to contract the disease, 
but more highly connected than average nodes. 

A well known result in network science is that more heterogeneous networks are less resistant to the invasion of diseases. 
However, in the self-organized networks studied here the opposite effect is observed. In comparison to random networks the 
self-organized networks are both more resistant to the disease and more heterogeneous in connectivity. 

The evolved networks gain their resistance to the disease from correlations between intrinsic parameters and the node connec¬ 
tivity. The increased resistance thus arises directly from the heterogeneity. A comparable effect is not possible in networks of 
identical nodes. While the specific bifurcation structure of the present model may depend on modelling assumptions, the basic 
interplay between the connectivity and the invasion threshold arises from the fundamental physics of the spreading process and is 
thus likely to be a generic feature that is observed across many models. It thus appears plausible that also in real world epidemics 
some anti-correlation between the true susceptibility and the effective degree of agents will be induced. By concentrating more 
of the remaining links of the more resilient individually the network will generally become both more heterogeneous and more 
resistant to the disease. 

It is possible that in real networks the self-organized heterogeneity is relatively minor compared to network heterogeneity 
that is unrelated to the epidemic in question. In that case the phenomenon described here would still occur but be of lesser 
importance. The extent to which self-organized heterogeneity plays a role in real world diseases will certainly depend on the 
disease in question, and can probably only be assessed through further empirical work. 
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FIG. 5: (Color online) Self-organized heterogeneity. Comparison of the degree distributions of the initial unadapted network used in the first 
set of simulations and the ‘adapted" self-organized network. The self-organized network is significantly more heterogeneous. However, at the 
same time it is more resistant to disease invasion. Parameters: /? = 0.032, t/> a = 0.65, ipb = 0.05, cu = 0.2, p = 0.002, N = 10 s , K = 10 6 . 


For epidemiology the results reported here imply, that network heterogeneity should not be considered in isolation. If there is 
an underlying heterogeneity in the susceptibility to the disease then a heterogeneous network may be more resistant to disease 
invasion than its homogeneous counterpart. Moreover, a vaccination campaign that targets the most highly connected nodes may 
end up vaccinating the wrong people as these nodes may also have the strongest natural resistance against the disease. 

Perhaps most importantly, our results suggest that the invasibility of the network is not governed by the heterogenity of the 
network alone, but by the heterogeneity of effective disease risk, which takes both node degree and susceptibility into account. 


Methods 

Moment-closure approximation. The time evolution of the proportion of nodes that are infected and of type u can then be 
written as 


^[Iu] = -MJ U ] + (14) 

v 

The two terms of this equation capture the recovery of infected nodes and the infection of susceptible nodes, respectively. In 
the second term we used the symbol [S U I V ] to denote the number of links between a node of state S and type u and a node of 
state I and type v, normalized with respect to the total number of nodes. This quantity therefore has the dimension of links per 
node. For determining these link densities, we can write additional evolution equations, which are in turn depend on the density 
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of triplet chains of nodes of given type and state [X U Y V Z W ], this yields 

S u s v ] = Klf i([s u i v ] + [i„s v ]) - Kl p J2^ s uS v U + Ms v Sui w }) 

( 15 ) 

+ ^£[S.][^4] + [S»][S V I W ]), 

2-^w w 


, \Iu ly ] - /,, /,;] + K% Pi'tpV -f- ti'y, ],S';, I v ] ) 

T ^ ^ [ /)< 5\i /,,; ] T 'lfj u ^ ] [/„ 5); I'll] ] ) ; 

it; it; 


(16) 


dt 


[^u-A;] —t?2/r[/y/y (/^ T [^u-A;] 

+ Pfa ^[S U S„/ W ] - ^ - u l s * I '> 


(17) 


where 


Kl 



U = V 

and «2 

u ^ v 


2 u = v 
1 u ^ v 


(18) 


In the equation for [S^SA], the S-S-links, the first term ki/^([jS u /„]+ [/„£>',,]) accounts for the creation of /y-SA-links by recov¬ 
ery of an infected node in an S-I-link. The factor k \ needs to be included to avoid double counting in the case of identical indices. 
The second term accounts for the destruction of .S'„-,S' r -Iinks due to infection of one of the two S-nodes by a node that is external 
to the link. The number of such infected nodes outside the S-S-link is given by the number of triplets [5' U S'„/ U) ] and [.S',,.S',; /,, ]. 
For a single S^-Sb-link the infection rate due to infected connected to the S^-node is /3i/jb{[S a SbIa\ + [<S a Sb7b])/[5'ai5'b]. he. the 
expected number of triplet chains that run through one specific the S-S-link and include an infected node on the S a side, mul¬ 
tiplied by the effective infection rate /3 1 /> 6 . To obtain the total rate we multiply by the number of those links, [S^S),], which 
cancels the denominator, take the infected on the other side of the S-S-link into account, and multiply k\ to avoid double¬ 
counting. The final term in the equation accounts for creation of S-S-links by rewiring. Links of the type S a -h, are rewired 
at rate to. When rewiring the susceptible node cuts the link to the infected and connects to a randomly chosen susceptible. In 
such a rewiring event, a new ,S' a -,S'„-linl< is created if the newly-chosen partner is of type A. This is the case with probability 
[£>]/([£„] + [Sk])• Such that the total rate of S^-Sa-link creation from S^-ib-link rewiring is uj[S a Ib] [<S' a ]/([S' a ] + [Sb]). Taking 
all possible combinations of indices and double-counting into account leads to the term in the equation. 

In the equation for /,, --1 i n ks the first term accounts for the loss of these links due to recovery of one of the linked nodes. 
The second term, k\ f3(ip v [I u S v \ + ipuiSuIv}), accounts for the creation of these links due to transmission of infection inside an 
S-I-link. Similarly, the final term accounts for the creation of I-I-link by infection of the S-node in an S-I-link from an internal 
source. In this term triplet chains appear in analogy to the term for the loss of S-S-links due to infection from sources external to 
the link, discussed above. In the equation for the S-I-links, the first term accounts for the creation of these links due to recovery 
of one infected node in an I-I-link. The second term accounts for the loss of S-I-links, both due to recovery of the I-node and 
due to internal transmission of the infection, resulting in an I-I-link. The fourth term captures the creation of S-I-links due to 
infection of an S-node in an S-S-link, whereas the fourth term captures the loss of S-I-links due to infection of the S-node from 
a source external to the link. Finally, the last term accounts for the loss of S-I-links due to rewiring. 

To cut the progression to ever larger network motifs one approximates the density of triplet chains by a moment closure 
approximation, here the pair approximation 


[X u Y v Z w ] = S [XuY ^ vZw] (19) 

[j-v\ 

where S is a factor arising from symmetries, such that 5 = 4 if X u = Y v = Z w , 5 = 2 if either X u = Y v ± Z w or 

X^Y„ = Z„,^6 = UfxjY„?Z„. 

Inherent in the moment-closure approximation is the assumption that long-ranged correlations vanish, such that the densities 
of motifs beyond the cut-off conform to statistical expectations. This assumption is the main source of inaccuracies in this 
type of approximation The approximation can still be used to identify phase transitions in the adaptive SIS model as 
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the correlations associated with these are still captured. However, the approximation performs poorly in fragmentation-type 
transitions, found for instance in the adaptive voter model lf53ll . 

The additional equations from the moment closure approximation that are used in the continuation are 


d[S h S h ] 

dt 


=p[S h I h ] - 2/3M [ShS tl [S ] hIa 


[5 b 5 b ][5 b 4] 

[^b] 


) + 


u[S h 


[5a] + [5b] 


([5 b 4] + [5b4]), 


( 20 ) 


d[S a S h } 
dt 


=/i([5 b / a ] + [5 a J b ]) - PM 


[5 b 5 a ][5 a /a 


[5 a 


[5 b 5 a ][5 a 4k a , , [5 a 5 b ][S' b / i 

-) - w]' 


[5a 


[5 a 5 b ][5 b /b] 

[5 b ] 


) + 


w[5 b 


[5 a ] + [5 b 


■([5 a / a ] + [5 a 4]) 


w[5 a 


[5 a ] + [5 b ] 


[5 b ] 

([Sb/a] + [5b 4]), 


( 21 ) 


d[5 a /a] 

dt 


=2/U[44] - (ft + P4 >a 
,[5 a Ja][S a Ja] 


- PM- 


-w)[S a / a ]+2^a( 
[5a4][5 a 4] 


[5 a 5 a ][5a4 

[5a] 


[5 a 


[5 a 




[5 a 5 a ] [5 a 4], 
[5a] ^ 


( 22 ) 


d[ShIh] = 2 /r[ 44 ] - (M + Wb + w)[Sb 4 ] + 2 pM- ShSh][Shh 


dt 


[5 b 


[5 b 5 b ][5 b 4] 

[5b] 


) 


ai , [5 b 4][5 b 4] , [5b4][5 b 4], 

- Pw b(-fTH- 1 -fTTl-), 


[5b] 


[5 b 


(23) 


d[5 a 4] _ (ji + + w ) [5a j b ] + p^ h (l SaSh][ShIa 


dt 


[5 b 


[5 a 5 b ][Sb4] 

[5b] 


_ ^, a( [5a4][5a4] + [5a4][5a4] )) 


[5 a 


[5 a 


(24) 


d[5 b / a ] =/i[4Jb] _ (/i + + W )[5 b /a] + PM- SaSh][SJa 


dt 


[5 a 


[S a Sb][Sa4] 

[5a] 


) 


_ ^, o( [Sb/a] [Sb/a] + [Sb/a][S b /b] )) 


[5 b 


[5 b 


(25) 


d[IJa] = - 2 /i[ 44 ] + ^a[S a /a] + PM^M^ + 


dt 


[5 a 


[5 a 


(26) 


d[,M = - 2,444] + ms, 4] + + MM,. 


dt 


[5 b 


[5 b 


(27) 
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Figure legends 

Figure 1: Bifurcation diagram of the adaptive heterogeneous SIS model. Show is the stationary level of disease prevalence 
I* as a function of infectivity 3. When the infectivitiy is decreased the endemic state vanishes in a saddle node bifurcation 
(‘Persistence threshold’). The disease free state can be invaded by the epidemic, if the infectivity surpasses a point where a 
transcritical bifurcation occurs (‘Invasion threshold’). Agent-based simulation (circles) and equation-based continuation (lines) 
provide consistent results on the persistence threshold, but predict different invasion thresholds. This discrepancy appears due 
to a projection effect, because the state where the disease is extinct is not uniquely defined (see text). In addition to stable 
solution branches (solid) the continuation also reveals an unstable solution branch (dotted). Parameters: 'd> a , = 0.65, ipb = 0.05, 
p a = 0.75, p = 0.002, N = 10 5 , K = 10 6 . 

Figure 2: Comparison of thresholds. The plot shows a very good agreement agreement between equation-based continuation 
(lines) and agent-based simulations (symbols) for the persistence thresholds /3 per (box, dashed). However, a notable difference 
exists for the invasion thresholds /3j nv (circle, dotted). Parameters: t/>b = 0.05, p = 0.002, N = 10 5 , K = 10 6 . 

Figure 3: Network heterogeneity in the adapted state, indicated by the degree ratio kb/k a . The dependence of the 
degree ratio in agent-based simulations (black circles) closely follows the prediction from integration of the ODE model (solid 
red line), a very good approximation is also provided by the relationship k\,/k a = \J'3 a /3)\ : , (blue dashed line), whereas the naive 
expectation kb/k a = V’a/V’b (pink dotted line) overestimates the network heterogeneity significantly. Contrary to expectations 
the networks following the naive solution (the most heterogeneous case), would be maximally stable against disease invasion. 
Parameters: if) b = 0.05, p = 0.002, N = 10, K = 10 6 . Inset: the magnification of the superimposed part of the main figure. 

Figure 4: Comparison of thresholds in self-organized networks. The plot shows a very good agreement between equation- 
based continuation (lines) and agent-based simulations started in an artificially created adapted state (symbols) for both the 
invasion thresholds /3j nv (circle, dotted) and the persistence thresholds /3 per (box, dashed). See Fig. 2 for comparison. Parameters: 
V> b = 0.05, p = 0.002, N = 10 5 , K = 10 6 . 

Figure 5: Self-organized heterogeneity. Comparison of the degree distributions of the initial unadapted network used in the 
first set of simulations and the ‘adapted’ self-organized network. The self-organized network is significantly more heterogeneous. 
However, at the same time it is more resistant to disease invasion. Parameters: /? = 0.032, ip a = 0.65, t/>b = 0.05, p = 0.002, 
N = 10 5 , K = 10 6 . 



